library(tidyverse)
library(patchwork)

dat1 = read.csv(file.choose()) # SSP245

dat2 = read.csv(file.choose()) # SSP585

dat3 = rbind(dat1,dat2)

# for Arabica remove elevation below 1000m.a.s.l
# for Robusta remove elevation above 14000m.a.s.l
dat3 = dat3 %>% filter(elevation>=1000)
dat3 = dat3 %>% filter(elevation<=1400)

current = dat3 %>% filter(scenario=='Current') %>% dplyr::select(-c(scenario))

p1 = ggplot(data=filter(dat3,scenario!='Current'),aes(elevation,suitability,color=period))+
  geom_smooth(data=current,aes(elevation,suitability))+
  geom_smooth()+
  facet_wrap(~scenario)+
  theme_bw()+
  labs(x='Elevation (m.a.s.l)',y='Suitability scores',title = '(A)')+
  theme(legend.position = 'none')

p2 = ggplot(data=filter(dat3,scenario!='Current'),aes(Y,suitability,color=period))+
  geom_smooth(data=current,aes(Y,suitability))+
  geom_smooth()+
  facet_wrap(~scenario)+
  theme_bw()+
  labs(x='Latitude',y='Suitability scores',title = '(B)',color='')+
  theme(legend.position = 'bottom')
  #scale_x_continuous(breaks = seq(-11.25,-10.25,0.3))# for Arabica

p3 = ggplot(data=filter(dat3,scenario!='Current'),aes(X,suitability,color=period))+
  geom_smooth(data=current,aes(X,suitability))+
  geom_smooth()+
  facet_wrap(~scenario)+
  theme_bw()+
  labs(x='Longitude',y='Suitability scores',title = '(C)')+
  theme(legend.position = 'none')

p4 = p1+p2+p3+plot_layout(ncol=2)

p4

ggsave('Robusta_lat_long_ele_suitability_change_ssp245_ssp585.png',p4,height=6,width=11,units='in',dpi=320)

rm(list=ls())
dev.off()
